Setup

library(easypackages)
libraries("here","ggplot2","reshape","ggpubr","ggrepel","forcats","patchwork","png","grid")

datapath = here("github_repo","data","recurrent_model")
figpath = here("github_repo","figures")

Figure 1

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

legend_labels_PSD = c(
                     "Piecewise fitting",
                      expression(g == 5.6),
                      expression(g == 14.8)
                     )

legend_labels_slope = c(expression(paste(nu[0] == 1.5," sp/s")),
                        expression(paste(nu[0] == 2," sp/s")),
                        'Baseline')

# Load data 
## LFPs
LFP_1 = read.csv(file.path(datapath,"LFP_5.646067415730337.csv"))
LFP_2 = read.csv(file.path(datapath,"LFP_14.79863985807809.csv"))

# PSDs
PSD_1 = read.csv(file.path(datapath,"PSD_5.646067415730337.csv"))
PSD_2 = read.csv(file.path(datapath,"PSD_14.79863985807809.csv"))
PSD_slope_1_lf = read.csv(file.path(datapath,"PSD_slope_1_5.646067415730337.csv"))
PSD_slope_1_hf = read.csv(file.path(datapath,"PSD_slope_2_5.646067415730337.csv"))
PSD_slope_2_lf = read.csv(file.path(datapath,"PSD_slope_1_14.79863985807809.csv"))
PSD_slope_2_hf = read.csv(file.path(datapath,"PSD_slope_2_14.79863985807809.csv"))

# LFP slopes 
lf_slopes_g_rate_1 = read.csv(file.path(datapath,"lf_slopes_g_rate_1.5.csv"))
lf_slopes_g_rate_2 = read.csv(file.path(datapath,"lf_slopes_g_rate_2.0.csv"))
hf_slopes_g_rate_1 = read.csv(file.path(datapath,"hf_slopes_g_rate_1.5.csv"))
hf_slopes_g_rate_2 = read.csv(file.path(datapath,"hf_slopes_g_rate_2.0.csv"))

# H values
LFP_H_g_rate_1 = read.csv(file.path(datapath,"LFP_H_g_rate_1.5.csv"))
LFP_H_g_rate_2 = read.csv(file.path(datapath,"LFP_H_g_rate_2.0.csv"))

# Baseline lines 
baseline_g_slopes_1 <- data.frame("var" = 11.3,"metric"=seq(-1, 1, length.out=100),"variable"='Baseline')
baseline_g_slopes_2 <- data.frame("var" = 11.3,"metric"=seq(-2.75, -1, length.out=100),"variable"='Baseline')
baseline_g_H <- data.frame("var" = 11.3,"metric"=seq(1.35,1.65, length.out=100),"variable"='Baseline')



# Figure 1A -------------------------------------------------------------------
recurrent_model_pdf = file.path(datapath,"recurrent_model.png")
img = readPNG(recurrent_model_pdf)
g = rasterGrob(img,interpolate = TRUE)

p = ggplot() + geom_blank()
p = p + annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) + 
  # theme_bw(panel.border = element_blank())
  theme_void()
p1a = p
p

# Figure 1B -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
LFP_1["variable"] <- "g=5.65" #'LFP 1'
LFP_2["variable"] <- "g=14.8" #'LFP 2'

LFP_all = rbind(LFP_1,LFP_2)
LFP_all$variable = factor(LFP_all$variable, levels = c("g=5.65","g=14.8"))

p1b = ggplot(data = LFP_all, aes(x = time, y = LFP, colour=variable)) + 
  facet_grid(variable ~ .) + 
  geom_line(size=0.25) + 
  scale_colour_manual(values = c("dodger blue","deep sky blue")) + guides(colour=FALSE) + 
  ylim(-2,10) + 
  ylab(" ") + 
  xlab("Time (ms)") + 
  ggtitle("LFP time series") + 
  theme(plot.title = element_text(hjust = 0.5, size=fontSize),
        strip.text.y = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize-2),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize))

# ggsave(filename = file.path(figpath,sprintf("Fig1B.pdf",fstem_oof)),plot=p1b)
p1b

# Figure 1C -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
PSD_1["variable"] <- 'PSD 1'
PSD_2["variable"] <- 'PSD 2'
PSD_slope_1_lf["variable"] <- 'PSD_slope_1_lf'
PSD_slope_1_hf["variable"] <- 'PSD_slope_1_hf'
PSD_slope_2_lf["variable"] <- 'PSD_slope_2_lf'
PSD_slope_2_hf["variable"] <- 'PSD_slope_2_hf'

merged_data_1 <- rbind(PSD_1,PSD_2)

p = ggplot()
p = p + geom_line(data = merged_data_1, aes(x = fx, y = PSD, colour = variable),size = 1) 
p = p + geom_line(data = PSD_slope_1_lf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_1_hf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_2_lf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_2_hf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 

p = p + xlab("Frequency (Hz)")+ ylab( expression(paste( mV^2/Hz)) ) + ggtitle("LFP PSD")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + theme(legend.title = element_blank(),legend.position = c(0.25, 0.25),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("dashed","solid","solid"),size=c(0.5,1,1))))
p = p + guides(colour = FALSE)
p = p + scale_colour_manual(values = c("black","dodger blue","deep sky blue"),labels = legend_labels_PSD)
p = p + scale_x_log10() + scale_y_log10()
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))
p1c = p
# ggsave(filename = file.path(figpath,sprintf("Fig1C.pdf",fstem_oof)))
p

# Figure 1D -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
lf_slopes_g_rate_1["variable"] <- '1.5 spikes/s'
lf_slopes_g_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(lf_slopes_g_rate_1,lf_slopes_g_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_g_slopes_1, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_g_slopes_1, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(20, 2, by = -2),1))
p = p + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("1/f slope") + ggtitle("Low-frequencies <30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + guides(colour=FALSE)
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

p1d = p
# ggsave(filename = file.path(figpath,sprintf("Fig1D.pdf",fstem_oof)))
p

# Figure 1E -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
hf_slopes_g_rate_1["variable"] <- '1.5 spikes/s'
hf_slopes_g_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(hf_slopes_g_rate_1,hf_slopes_g_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_g_slopes_2, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_g_slopes_2, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(20, 2, by = -2),1))
p = p + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("1/f slope") + ggtitle("High-frequencies >30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p1e = p
# ggsave(filename = file.path(figpath,sprintf("Fig1E.pdf",fstem_oof)))
p

# Figure 1F -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
LFP_H_g_rate_1["variable"] <- '1.5 spikes/s'
LFP_H_g_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(LFP_H_g_rate_1,LFP_H_g_rate_2)

p1 = ggplot()
p1 = p1 + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p1 = p1 + geom_line(data = baseline_g_H, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_g_H, aes(x = var, y = metric, colour = FALSE))

p1 = p1 + scale_x_reverse(breaks = round(seq(20, 2, by = -2),1))
p1 = p1 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ggtitle("Hurst exponent (H)")
p1 = p1 + theme(plot.title = element_text(hjust = 0.5))
# p1 = p1 + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p1 = p1 + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p1 = p1 + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p1 = p1 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")
# ggsave(filename = file.path(figpath,sprintf("Fig1Fa.pdf",fstem_oof)), plot=p1)

# Box plot
box1 <- data.frame(seq(1,60,1),seq(1,60,1))
box2 <- data.frame(seq(1,60,1),seq(1,60,1))

# Clustering of H in 3 groups
box1[1:20,2] = LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] < 7.5,3]
box1[21:40,2] = LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] >= 7.5 & LFP_H_g_rate_1[,"var"] < 10,3]
box1[41:60,2] = LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] >= 11,3]
box1[1:20,1] = format(round(mean(LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] < 7.5,2]), 2), nsmall = 2)
box1[21:40,1] = format(round(mean(LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] >= 7.5 & LFP_H_g_rate_1[,"var"] < 10,2]), 2), nsmall = 2)
box1[41:60,1] = format(round(mean(LFP_H_g_rate_1[LFP_H_g_rate_1[,"var"] >= 11,2]), 2), nsmall = 2)
colnames(box1) <- c("var","metric")

box2[1:20,2] = LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 5.8 & LFP_H_g_rate_2[,"var"] < 7.5,3]
box2[21:40,2] = LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 7.8 & LFP_H_g_rate_2[,"var"] < 10,3]
box2[41:60,2] = LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 11 & LFP_H_g_rate_2[,"var"] < 14.6,3]
box2[1:20,1] = format(round(mean(LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 5.8 & LFP_H_g_rate_2[,"var"] < 7.5,2]), 2), nsmall = 2)
box2[21:40,1] = format(round(mean(LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 7.8 & LFP_H_g_rate_2[,"var"] < 10,2]), 2), nsmall = 2)
box2[41:60,1] = format(round(mean(LFP_H_g_rate_2[LFP_H_g_rate_2[,"var"] >= 11 & LFP_H_g_rate_2[,"var"] < 14.6,2]), 2), nsmall = 2)
colnames(box2) <- c("var","metric")

# boxplot
p2 <-ggplot(box1, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p2 = p2 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p2 = p2 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p2 = p2 + annotate(geom="text", x=3, y=1.6, label="1.5 sp/s", color="red", size = 5)

# Other props.
p2 = p2 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.4,1.65) 
#p2 = p2 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
# p2 = p2 + scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9")) + ggtitle("Hurst exponent (H)")
p2 = p2 + scale_colour_manual(values = c("#808000","#808000","#808000")) + ggtitle("1.5 sp/s")
p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust = 0.5),
        legend.position = "none")
# ggsave(filename = file.path(figpath,sprintf("Fig1Fb.pdf",fstem_oof)), plot=p2)

# boxplot
p3 <-ggplot(box2, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p3 = p3 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p3 = p3 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p3 = p3 + annotate(geom="text", x=3, y=1.5, label="2 sp/s", color="red", size = 5)

# Other props.
p3 = p3 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.4,1.55)
#p3 = p3 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
# p3 = p3 + scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9")) + ggtitle("Hurst exponent (H)")
p3 = p3 + scale_colour_manual(values = c("#FA8072","#FA8072","#FA8072")) + ggtitle("2 sp/s")
p3 = p3 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust = 0.5),
        legend.position = "none")
# ggsave(filename = file.path(figpath,sprintf("Fig1Fc.pdf",fstem_oof)), plot=p3)

p1f = p1 / (p2 | p3)
# multiplot(p1, p2, p3,layout = matrix(c(1,1,2,3),nrow=2, byrow=TRUE))
# ggsave(filename = file.path(figpath,sprintf("Fig1F.pdf",fstem_oof)), plot=p1f)
p1f

# Put plots together
p_final = (p1a | p1b | p1c) / (p1d | p1e | p1) / (p2 | p3) + plot_layout(heights =c(1,1,1.1))
ggsave(filename = file.path(figpath,sprintf("Fig1.pdf",fstem_oof)), plot=p_final)
p_final

Figure 2

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

# Load data 
## HRF
HRF = read.csv(file.path(datapath,"HRF.csv"))

## FFTs of HRF and HPF
fft_HRF = read.csv(file.path(datapath,"fft_HRF.csv"))
fft_HPF = read.csv(file.path(datapath,"fft_HPF.csv"))

# H values
BOLD_H_g_rate_1 = read.csv(file.path(datapath,"BOLD_H_g_rate_1.5_Magri_kernel_Yes_HPF.csv"))
BOLD_H_g_rate_2 = read.csv(file.path(datapath,"BOLD_H_g_rate_2.0_Magri_kernel_Yes_HPF.csv"))

# Correlation of LFP power and BOLD 
corr_matrix = read.csv(file.path(datapath,"corr_matrix_Magri_kernel_Yes_HPF.csv"))
alpha_band = read.csv(file.path(datapath,"alpha_band_Magri_kernel_Yes_HPF.csv"))
beta_band = read.csv(file.path(datapath,"beta_band_Magri_kernel_Yes_HPF.csv"))
gamma_band = read.csv(file.path(datapath,"gamma_band_Magri_kernel_Yes_HPF.csv"))
LFP_band = read.csv(file.path(datapath,"LFP_band_Magri_kernel_Yes_HPF.csv"))

# Baseline lines 
baseline_g_H <- data.frame("var" = 11.3,"metric"=seq(1.2,1.65, length.out=100),"variable"='Baseline')


# Figure 2A -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

colnames(corr_matrix) <- c('X.1','X','Y','r')

# jet colormap
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

cMid = 0.1
cLims = c(-0.1,0.3)

p = ggplot(corr_matrix, aes(X, Y, fill= r)) + geom_tile()
# p = p + scale_fill_gradientn(colors = jet.colors(7))
# p = p + scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100))
p = p + scale_fill_gradient2(low="blue", mid="white", high="red",
                             midpoint=cMid, limits=cLims)
p = p + xlab("Delay (s)") + ylab("Frequency (Hz)")+ ggtitle("Correlation between LFP power and BOLD")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + scale_y_reverse(breaks = round(seq(100, 0, by = -20),1))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.title = element_text(size=fontSize,hjust=0.25))

p2a = p
# ggsave(filename = file.path(figpath,sprintf("Fig2A.pdf",fstem_oof)))
p

# Figure 2B -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
alpha_band["variable"] <- 'alpha'
beta_band["variable"] <- 'beta'
gamma_band["variable"] <- 'gamma'
LFP_band["variable"] <- 'LFP'

merged_data <- rbind(alpha_band,beta_band,gamma_band,LFP_band)

p = ggplot()
p = p + geom_line(data = merged_data, aes(x = time, y = band, colour = variable),size = 1) 
p = p + xlab("Delay (s)") + ylab("Correlation") + ggtitle("Correlation in selected frequency bands")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + theme(legend.title = element_blank(),legend.position = c(0.9, 0.85),legend.text = element_text(size=fontSize),legend.text.align = 0)
p = p + scale_colour_manual(values = c("#006400","red","blue","black")) + guides(colour=FALSE)

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

p2b = p
# ggsave(filename = file.path(figpath,sprintf("Fig2B.pdf",fstem_oof)))
p

# Figure 2C and 2D ------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

# Box plot
box1 <- data.frame(seq(1,60,1),seq(1,60,1))
box2 <- data.frame(seq(1,60,1),seq(1,60,1))

# Clustering of H in 3 groups
box1[1:20,2] = BOLD_H_g_rate_1[seq(0,60,3),3]
box1[21:40,2] = BOLD_H_g_rate_1[seq(1,60,3),3]
box1[41:60,2] = BOLD_H_g_rate_1[seq(2,60,3),3]
box1[1:20,1] = format(round(mean(BOLD_H_g_rate_1[seq(0,60,3),2]), 2), nsmall = 2)
box1[21:40,1] = format(round(mean(BOLD_H_g_rate_1[seq(1,60,3),2]), 2), nsmall = 2)
box1[41:60,1] = format(round(mean(BOLD_H_g_rate_1[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box1) <- c("var","metric")

box2[1:20,2] = BOLD_H_g_rate_2[seq(0,60,3),3]
box2[21:40,2] = BOLD_H_g_rate_2[seq(1,60,3),3]
box2[41:60,2] = BOLD_H_g_rate_2[seq(2,60,3),3]
box2[1:20,1] = format(round(mean(BOLD_H_g_rate_2[seq(0,60,3),2]), 2), nsmall = 2)
box2[21:40,1] = format(round(mean(BOLD_H_g_rate_2[seq(1,60,3),2]), 2), nsmall = 2)
box2[41:60,1] = format(round(mean(BOLD_H_g_rate_2[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box2) <- c("var","metric")

# boxplot
p2 <-ggplot(box1, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p2 = p2 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p2 = p2 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p2 = p2 + annotate(geom="text", x=3, y=1.6, label="1.5 sp/s", color="red", size = 6)

# Other props.
p2 = p2 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.2,1.65) + ggtitle("1.5 sp/s")
p2 = p2 + theme(plot.title = element_text(hjust = 0.5))
#p2 = p2 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p2 = p2 + scale_color_manual(values=c("#808000", "#808000", "#808000"))
p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p2c = p2
# ggsave(filename = file.path(figpath,sprintf("Fig2C.pdf",fstem_oof)),plot=p2c)
p2

# boxplot
p3 <-ggplot(box2, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p3 = p3 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p3 = p3 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p3 = p3 + annotate(geom="text", x=3, y=1.5, label="2 sp/s", color="red", size = 6)

# Other props.
p3 = p3 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.1,1.55) + ggtitle("2 sp/s")
p3 = p3 + theme(plot.title = element_text(hjust = 0.5))
#p3 = p3 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p3 = p3 + scale_color_manual(values=c("#FA8072", "#FA8072", "#FA8072"))
p3 = p3 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p2d = p3
# ggsave(filename = file.path(figpath,sprintf("Fig2D.pdf",fstem_oof)),plot=p2d)
p3

# Put plots together
p_final = (p2a | p2b) / (p2c | p2d)
ggsave(filename = file.path(figpath,sprintf("Fig2.pdf",fstem_oof)), plot=p_final)
p_final

Figure 3

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

legend_labels_slope = c(expression(paste(nu[0] == 1.5," sp/s")),
                        expression(paste(nu[0] == 2," sp/s")),
                        'Baseline')


# Data for Figure 3A-B
# H values
LFP_H_vth_rate_1 = read.csv(file.path(datapath,"LFP_H_vth_rate_1.5.csv"))
LFP_H_vth_rate_2 = read.csv(file.path(datapath,"LFP_H_vth_rate_2.0.csv"))

# H values
LFP_H_EL_rate_1 = read.csv(file.path(datapath,"LFP_H_EL_rate_1.5.csv"))
LFP_H_EL_rate_2 = read.csv(file.path(datapath,"LFP_H_EL_rate_2.0.csv"))

# Baseline lines 
baseline_vth_H <- data.frame("var" = -52,"metric"=seq(1.35,1.7, length.out=100),"variable"='Baseline')

baseline_EL_H <- data.frame("var" = -70,"metric"=seq(1.35,1.75, length.out=100),"variable"='Baseline')


# Figure 3A -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
LFP_H_vth_rate_1["variable"] <- '1.5 spikes/s'
LFP_H_vth_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(LFP_H_vth_rate_1,LFP_H_vth_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_vth_H, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_vth_H, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( V[th]," (mV)")) ) + ylab("H") + ggtitle("Hurst exponent (H)")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072","black"),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p3a = p
# ggsave(filename = file.path(figpath,sprintf("Fig3A.pdf",fstem_oof)))
p

# Figure 3B -------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
LFP_H_EL_rate_1["variable"] <- '1.5 spikes/s'
LFP_H_EL_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(LFP_H_EL_rate_1,LFP_H_EL_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_EL_H, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_EL_H, aes(x = var, y = metric, colour = FALSE))
# p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( E[L]," (mV)")) ) + ylab("H") + ggtitle("Hurst exponent (H)")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p3b = p
# ggsave(filename = file.path(figpath,sprintf("Fig3B.pdf",fstem_oof)))
p

# Figure 3C -------------------------------------------------------------------
# initial information
window_size = 512
nwindows = 3529
vol_inject_start = 970
vol_treatment_start = vol_inject_start+900 

# read in sliding-window analysis data
datapath_dreadd = here("github_repo","data","dreadd")
data_win = read.csv(file.path(datapath_dreadd,"pheno_data+Hwin_dreaddexcitation.csv"))

# melt the data frame into long format
data_win_long = melt(data_win, 
                     id.vars = c("filename","condition","scan_day","H"))
colnames(data_win_long)[ncol(data_win_long)] = "Hwin"

# make a time column
data_win_long$time = NA
for (i in 1:nwindows){
  mask = data_win_long$variable==sprintf("window_%04d",i)
  data_win_long[mask,"time"] = i
}

# make column indicating treatment phase
data_win_long$treatment_phase = NA
data_win_long$treatment_phase[data_win_long$time<=(vol_inject_start-window_size)-1] = "Baseline"
data_win_long$treatment_phase[data_win_long$time>=(vol_inject_start-window_size) & data_win_long$time<=(vol_treatment_start-window_size)-1] = "Transition"
data_win_long$treatment_phase[data_win_long$time>(vol_treatment_start-window_size)] = "Treatment"
data_win_long$treatment_phase = factor(data_win_long$treatment_phase)


# plotting using GAM smoother for both individual and group trajectories
yLims = c(1.4,1.8)
p = ggplot(data = data_win_long, 
           aes(x = time,
               y = Hwin,
               group = filename)) +
            facet_grid(. ~ condition)

# plot data, one line per mouse
# p = p + geom_line(alpha=0.5,colour="gray75")

# add smoothed individual lines for each mouse
p = p + geom_smooth(se = FALSE ,colour='gray75',alpha=0.1)

# add group-level fit line
p = p + geom_smooth(se=TRUE, aes(group = interaction(condition,treatment_phase), colour = treatment_phase), size=4)

# add other stuff to the plot
p = p + geom_vline(xintercept = (vol_inject_start-window_size)) +
  geom_vline(xintercept = vol_treatment_start-window_size) +
  ylab("H") +
  xlab("Time Window") + ggtitle("DREADD Excitation") +
  scale_y_continuous(limits = yLims) + guides(colour=FALSE)

p = p + theme(axis.text.x = element_text(size=fontSize),
              axis.text.y = element_text(size=fontSize),
              axis.title.x = element_text(size=fontSize),
              axis.title.y = element_text(size=fontSize),
              strip.text.x = element_text(size=fontSize),
              plot.title = element_text(size=fontSize,hjust = 0.5))

p3c = p
# save figure
# ggsave(filename = file.path(figpath,"Fig3C.pdf"))
p

# Figure 3D -------------------------------------------------------------------
# initial information
window_size = 512
nwindows = 3389
vol_inject_start = 900
vol_treatment_start = 1800

# read in sliding-window analysis data
data_win = read.csv(file.path(datapath_dreadd,"pheno_data+Hwin_dreaddsilencing.csv"))

# melt the data frame into long format
data_win_long = melt(data_win, 
                     id.vars = c("filename","condition","scan_day","H"))
colnames(data_win_long)[ncol(data_win_long)] = "Hwin"

# make a time column
data_win_long$time = NA
for (i in 1:nwindows){
  mask = data_win_long$variable==sprintf("window_%04d",i)
  data_win_long[mask,"time"] = i
}

# make column indicating treatment phase
data_win_long$treatment_phase = NA
data_win_long$treatment_phase[data_win_long$time<=(vol_inject_start-window_size)-1] = "Baseline"
data_win_long$treatment_phase[data_win_long$time>=(vol_inject_start-window_size) & data_win_long$time<=(vol_treatment_start-window_size)-1] = "Transition"
data_win_long$treatment_phase[data_win_long$time>(vol_treatment_start-window_size)] = "Treatment"
data_win_long$treatment_phase = factor(data_win_long$treatment_phase)

# plotting using GAM smoother for both individual and group trajectories
yLims = c(1.4,1.8)
p = ggplot(data = data_win_long, 
           aes(x = time,
               y = Hwin,
               group = filename)) +
            facet_grid(. ~ condition)

# plot data, one line per mouse
# p = p + geom_line(alpha=0.5,colour="gray75")

# add smoothed individual lines for each mouse
p = p + geom_smooth(se = FALSE ,colour='gray75',alpha=0.1)

# add group-level fit line
p = p + geom_smooth(se=TRUE, aes(group = interaction(condition,treatment_phase), colour = treatment_phase), size=4)

# add other stuff to the plot
p = p + geom_vline(xintercept = (vol_inject_start-window_size)) +
  geom_vline(xintercept = vol_treatment_start-window_size) +
  ylab("Hurst Exponent (H)") +
  xlab("Time Window") + ggtitle("DREADD Silencing") + 
  scale_y_continuous(limits = yLims) + guides(colour=FALSE)

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust = 0.5))

p3d = p
# save figure
# ggsave(filename = file.path(figpath,"Fig3D.pdf"))
p

# Put plots together
p_final = (p3a | p3b) / (p3c | p3d)
ggsave(filename = file.path(figpath,sprintf("Fig3.pdf",fstem_oof)), plot=p_final)
p_final

Supp Figure 1

fontSize = 8
datapath_gao = here("github_repo","data","gao_model")

fstem_oof = "oof"
file2use = file.path(datapath_gao,sprintf("Hsim_%s.csv",fstem_oof))
data_oof = read.csv(file2use)

# H of BOLD
BOLD_H = read.csv(file.path(datapath,"Gao_BOLD_H_Magri_HRF.csv"))
data_oof[,4] = BOLD_H[,3]

# Supp Fig 1A -----------------------------------------------------------------
gao_model_pdf = file.path(datapath_gao,"gao_model.png")
img = readPNG(gao_model_pdf)
g = rasterGrob(img,interpolate = TRUE)

p = ggplot() + geom_blank()
p = p + annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) + 
  theme_void()
ps1a = p
p

# Supp Fig 1B -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

ei_oof_data = read.csv(file.path(datapath_gao,"EIsim_oof.csv"))
nchannels = 128
ei_oof_plot = ggplot(data = ei_oof_data, aes(x = EIratio, y = channel_001))
for (ichan in 1:nchannels){
  chan_name = sprintf("channel_%03d",ichan)
  ei_oof_plot = ei_oof_plot + geom_line(aes_string(y=chan_name), alpha = 0.3, colour = "grey60")  
}
ei_oof_data$mean_H = rowMeans(ei_oof_data[,3:ncol(ei_oof_data)])
ei_oof_plot = ei_oof_plot + geom_smooth(data = ei_oof_data, aes(y=ei_oof_data$mean_H))
ei_oof_plot = ei_oof_plot + scale_x_continuous(breaks=ei_oof_data$EIratio[labelidx2use],
                                               labels=ei_oof_data$EIratio[labelidx2use])
ei_oof_plot = ei_oof_plot + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("1/f slope") + ggtitle("Ratio of synaptic conductances")
ei_oof_plot = ei_oof_plot + theme(plot.title = element_text(hjust = 0.5))
ei_oof_plot = ei_oof_plot + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))
ei_oof_plot = ei_oof_plot + scale_x_reverse()

ps1b = ei_oof_plot
# ggsave(filename = file.path(figpath,"SuppFig1B.pdf"), plot=ei_oof_plot)
ei_oof_plot

# Supp Fig 1C -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

ei_data = read.csv(file.path(datapath_gao,"EIsim_H.csv"))
nchannels = 128
ei_H_plot = ggplot(data = ei_data, aes(x = EIratio, y = channel_001))
for (ichan in 1:nchannels){
  chan_name = sprintf("channel_%03d",ichan)
  ei_H_plot = ei_H_plot + geom_line(aes_string(y=chan_name), alpha = 0.3, colour = "grey60")  
}
ei_data$mean_H = rowMeans(ei_data[,3:ncol(ei_data)])
ei_H_plot = ei_H_plot + geom_smooth(data = ei_data, aes(y=ei_data$mean_H))
ei_H_plot = ei_H_plot + scale_x_continuous(breaks=ei_data$EIratio[labelidx2use],
                                           labels=ei_data$EIratio[labelidx2use])
ei_H_plot = ei_H_plot + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ggtitle("Ratio of synaptic conductances")
ei_H_plot = ei_H_plot + theme(plot.title = element_text(hjust = 0.5))
ei_H_plot = ei_H_plot + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))
ei_H_plot = ei_H_plot + scale_x_reverse()


ps1c = ei_H_plot
# ggsave(filename = file.path(figpath,"SuppFig1C.pdf"), plot=ei_H_plot)
ei_H_plot

# Supp Fig 1D -----------------------------------------------------------------
y_limits = c(1.4,2)
labelidx2use = c(1,6,11,16,21)

melted_data_oof = melt(data_oof, id.vars = c("legend_labels","oof"))

p = ggplot(data = melted_data_oof, aes(x = legend_labels, y = value, colour = variable))
p = p + geom_point() + geom_smooth() + guides(colour = FALSE)
p = p + scale_x_continuous(breaks=data_oof$legend_labels[labelidx2use], labels=data_oof$legend_labels[labelidx2use])
p = p + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
p = p + xlab("1/f slope in LFP data") + ylab("H") + ggtitle("Simulated LFP (1/f only)")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps1d = p
# ggsave(filename = file.path(figpath,"SuppFig1D.pdf"), plot=p)
p

# Put plots together
p_final = (ps1a | ps1b) / (ps1c | ps1d)
ggsave(filename = file.path(figpath,"SuppFig1.pdf"), plot=p_final)
p_final

Supp Figure 2

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

legend_labels_PSD = c(
                     "FOOOF fitting",
                      expression(g == 5.6),
                      expression(g == 14.8)
                     )

legend_labels_slope = c(expression(paste(nu[0] == 1.5," sp/s")),
                        expression(paste(nu[0] == 2," sp/s")),
                        'Baseline')

# Load data 
# PSDs
PSD_1 = read.csv(file.path(datapath,"PSD_2_5.646067415730337.csv"))
PSD_2 = read.csv(file.path(datapath,"PSD_2_14.79863985807809.csv"))
fooof_1 = read.csv(file.path(datapath,"fooof_5.646067415730337.csv"))
fooof_2 = read.csv(file.path(datapath,"fooof_14.79863985807809.csv"))

# Ap. slopes
fooof_ap_g_rate_1 = read.csv(file.path(datapath,"fooof_ap_g_rate_1.5.csv"))
fooof_ap_g_rate_2 = read.csv(file.path(datapath,"fooof_ap_g_rate_2.0.csv"))

# Baseline lines 
baseline_g_ap_slopes <- data.frame("var" = 11.3,"metric"=seq(-1.25, 0, length.out=100),"variable"='Baseline')

# Supp Fig 2A -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
PSD_1["variable"] <- 'PSD 1'
PSD_2["variable"] <- 'PSD 2'
fooof_1["variable"] <- 'PSD_slope_1_ap'
fooof_2["variable"] <- 'PSD_slope_2_ap'

merged_data_1 <- rbind(PSD_1,PSD_2)

p = ggplot()
p = p + geom_line(data = merged_data_1, aes(x = fx, y = PSD, colour = variable),size = 1) 
p = p + geom_line(data = fooof_1, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = fooof_2, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 

p = p + xlab("Frequency (Hz)")+ ylab( expression(paste( mV^2/Hz)) ) + ggtitle("LFP PSD")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + theme(legend.title = element_blank(),legend.position = c(0.25, 0.25),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("dashed","solid","solid"),size=c(0.5,1,1))))
p = p + guides(colour = FALSE)
p = p + scale_colour_manual(values = c("black","dodger blue","deep sky blue"),labels = legend_labels_PSD)
p = p + scale_x_log10() + scale_y_log10()
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps2a = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig2A.pdf",fstem_oof)))
p

# Supp Fig 2B -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
fooof_ap_g_rate_1["variable"] <- '1.5 spikes/s'
fooof_ap_g_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(fooof_ap_g_rate_1,fooof_ap_g_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_g_ap_slopes, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_g_ap_slopes, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(20, 2, by = -2),1))
p = p + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("1/f slope") + ggtitle("Aperiodic slope")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps2b = p 
# ggsave(filename = file.path(figpath,sprintf("SuppFig2B.pdf",fstem_oof)))
p

# Put plots together
p_final = (ps2a | ps2b)
ggsave(filename = file.path(figpath,sprintf("SuppFig2.pdf",fstem_oof)), plot=p_final)
p_final

Supp Figure 3

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

# Load data
# LFPs
LFP_1 = read.csv(file.path(datapath,"LFP_5.646067415730337.csv"))
LFP_2 = read.csv(file.path(datapath,"LFP_14.79863985807809.csv"))

# PSDs
PSD_1 = read.csv(file.path(datapath,"PSD_5.646067415730337.csv"))
PSD_2 = read.csv(file.path(datapath,"PSD_14.79863985807809.csv"))
PSD_slope_1_lf = read.csv(file.path(datapath,"PSD_slope_1_5.646067415730337.csv"))
PSD_slope_1_hf = read.csv(file.path(datapath,"PSD_slope_2_5.646067415730337.csv"))
PSD_slope_2_lf = read.csv(file.path(datapath,"PSD_slope_1_14.79863985807809.csv"))
PSD_slope_2_hf = read.csv(file.path(datapath,"PSD_slope_2_14.79863985807809.csv"))

## HRF
HRF = read.csv(file.path(datapath,"HRF.csv"))

## FFTs of HRF and HPF
fft_HRF = read.csv(file.path(datapath,"fft_HRF.csv"))
fft_HPF = read.csv(file.path(datapath,"fft_HPF.csv"))

## FFTs
fft_BOLD_1_HPF_No = read.csv(file.path(datapath,"fft_BOLD_HPF_No_6.333595920236343.csv"))
fft_BOLD_2_HPF_No = read.csv(file.path(datapath,"fft_BOLD_HPF_No_13.409410112357872.csv"))
fft_BOLD_1_HPF_Yes = read.csv(file.path(datapath,"fft_BOLD_HPF_Yes_6.333595920236343.csv"))
fft_BOLD_2_HPF_Yes = read.csv(file.path(datapath,"fft_BOLD_HPF_Yes_13.409410112357872.csv"))

# BOLD
BOLD_1 = read.csv(file.path(datapath,"BOLD_6.333595920236343.csv"))
BOLD_2 = read.csv(file.path(datapath,"BOLD_13.409410112357872.csv"))

# Supp Fig 3A -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
LFP_1["variable"] <- "g=5.65" #'LFP 1'
LFP_2["variable"] <- "g=14.8" #'LFP 2'

LFP_all = rbind(LFP_1,LFP_2)
LFP_all$variable = factor(LFP_all$variable, levels = c("g=5.65","g=14.8"))

p1b = ggplot(data = LFP_all, aes(x = time, y = LFP, colour=variable)) + 
  facet_grid(variable ~ .) + 
  geom_line(size=0.5) + 
  scale_colour_manual(values = c("dodger blue","deep sky blue")) + guides(colour=FALSE) + 
  ylim(-2,10) + 
  ylab(" ") + 
  xlab("Time (ms)") + 
  ggtitle("LFP time series") + 
  theme(plot.title = element_text(hjust = 0.5, size=fontSize),
        strip.text.y = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize-2),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize))

ps3a = p1b
# ggsave(filename = file.path(figpath,sprintf("SuppFig3A.pdf",fstem_oof)),plot=p1b)
p1b

# Supp Fig 3B -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
HRF["variable"] <- 'HRF'
fft_HRF["variable"] <- 'fft_HRF'
fft_HPF["variable"] <- 'fft_HPF'

merged_data_2 <- rbind(fft_HRF,fft_HPF)

p1 = ggplot()
p2 = ggplot()

p1 = p1 + geom_line(data = HRF, aes(x = time, y = HRF, colour = 'black'),size = 1) 
p2 = p2 + geom_line(data = merged_data_2, aes(x = fx, y = PSD, colour = variable),size = 1) 

p1 = p1 + xlab("Time (ms)") + ggtitle("HRF")
p1 = p1 + theme(plot.title = element_text(hjust = 0.5)) + guides(color=FALSE)

p2 = p2 + xlab("Frequency (Hz)") + ggtitle("FFTs of convolution kernels")
p2 = p2 + theme(plot.title = element_text(hjust = 0.5)) + guides(color=FALSE)

# p1 = p1 + theme(legend.title = element_blank(),legend.position = c(0.85, 0.4),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p2 = p2 + theme(legend.title = element_blank(),legend.position = c(0.85, 0.4),legend.text = element_text(size=fontSize),legend.text.align = 0)

p1 = p1 + scale_colour_manual(values = c("black"))
p2 = p2 + scale_colour_manual(values = c("red","blue"))#, labels = c("HPF","HRF"))

p2 = p2 + scale_x_log10() + scale_y_log10()

p1 = p1 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_blank(),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_blank(),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps3b = (p1 / p2)
# ggsave(filename = file.path(figpath,sprintf("SuppFig3B.pdf",fstem_oof)),ps3b)
ps3b

# Supp Fig 3C -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
PSD_1["variable"] <- 'PSD 1'
PSD_2["variable"] <- 'PSD 2'
PSD_slope_1_lf["variable"] <- 'PSD_slope_1_lf'
PSD_slope_1_hf["variable"] <- 'PSD_slope_1_hf'
PSD_slope_2_lf["variable"] <- 'PSD_slope_2_lf'
PSD_slope_2_hf["variable"] <- 'PSD_slope_2_hf'

merged_data_1 <- rbind(PSD_1,PSD_2)

p = ggplot()
p = p + geom_line(data = merged_data_1, aes(x = fx, y = PSD, colour = variable),size = 1) 
p = p + geom_line(data = PSD_slope_1_lf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_1_hf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_2_lf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 
p = p + geom_line(data = PSD_slope_2_hf, aes(x = fx, y = PSD, colour = "black"),size = 0.5,linetype = "dashed") 

p = p + xlab("Frequency (Hz)")+ ylab( expression(paste( mV^2/Hz)) ) + ggtitle("LFP PSD")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + theme(legend.title = element_blank(),legend.position = c(0.25, 0.25),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("dashed","solid","solid"),size=c(0.5,1,1))))
p = p + guides(colour = FALSE)
p = p + scale_colour_manual(values = c("black","dodger blue","deep sky blue"))#,labels = legend_labels_PSD)
p = p + scale_x_log10() + scale_y_log10()
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))
ps3c = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig3C.pdf",fstem_oof)))
p

# Supp Fig 3D -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
BOLD_1["variable"] <- 'BOLD_1'
BOLD_2["variable"] <- 'BOLD_2'

p1 = ggplot()
p2 = ggplot()

p1 = p1 + geom_line(data = BOLD_1, aes(x = time, y = BOLD, colour = variable),size = 1)
p2 = p2 + geom_line(data = BOLD_2, aes(x = time, y = BOLD, colour = variable),size = 1)

p1 = p1  + xlab("Time (s)") + ggtitle("Recovered BOLD signals") + ylim(-2,1) + guides(colour=FALSE)
p1 = p1 + theme(plot.title = element_text(hjust = 0.5))

p2 = p2 + xlab("Time (s)") + ylim(-2,1) + guides(colour=FALSE)
p2 = p2 + theme(plot.title = element_text(hjust = 0.5))

# p1 = p1 + theme(legend.title = element_blank(),legend.position = c(0.85, 0.25),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p2 = p2 + theme(legend.title = element_blank(),legend.position = c(0.85, 0.25),legend.text = element_text(size=fontSize),legend.text.align = 0)

p1 = p1 + scale_colour_manual(values = c("dodger blue"),labels = c(expression(g == 6.3)))
p2 = p2 + scale_colour_manual(values = c("deep sky blue"),labels = c(expression(g == 13.4)))

p1 = p1 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_blank(),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_blank(),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_blank())

ps3d = (p1 / p2)
# ggsave(filename = file.path(figpath,sprintf("SuppFig3D.pdf",fstem_oof)),ps3d)
ps3d

# Supp Fig 3E -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
fft_BOLD_1_HPF_No["variable"] <- 'fft_BOLD_1_HPF_No'
fft_BOLD_2_HPF_No["variable"] <- 'fft_BOLD_2_HPF_No'

merged_data <- rbind(fft_BOLD_1_HPF_No,fft_BOLD_2_HPF_No)

p1 = ggplot()

p1 = p1 + geom_line(data = merged_data, aes(x = fx, y = PSD, colour = variable),size = 1)
p1 = p1 + xlab("Frequency (Hz)")
p1 = p1 + ylab("Normalized FFTs of BOLD signals")
p1 = p1 + ggtitle("Without high-pass filter")
p1 = p1 + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p1 = p1 + theme(legend.title = element_blank(),legend.position = c(0.75, 0.65),legend.text = element_text(size=fontSize),legend.text.align = 0)
p1 = p1 + scale_colour_manual(values = c("dodger blue","deep sky blue"),labels = c(expression(g == 6.3),expression(g == 13.4)))

p1 = p1 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps3e = p1
# ggsave(filename = file.path(figpath,sprintf("SuppFig3E.pdf",fstem_oof)),plot=ps3d)
p1

# Supp Fig 3F -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
fft_BOLD_1_HPF_Yes["variable"] <- 'fft_BOLD_1_HPF_Yes'
fft_BOLD_2_HPF_Yes["variable"] <- 'fft_BOLD_2_HPF_Yes'

merged_data <- rbind(fft_BOLD_1_HPF_Yes,fft_BOLD_2_HPF_Yes)

p1 = ggplot()

p1 = p1 + geom_line(data = merged_data, aes(x = fx, y = PSD, colour = variable),size = 1)
p1 = p1 + xlab("Frequency (Hz)") + ylab("Normalized FFTs of BOLD signals")
p1 = p1 + ggtitle("With high-pass filter")
p1 = p1 + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p1 = p1 + theme(legend.title = element_blank(),legend.position = c(0.75, 0.65),legend.text = element_text(size=fontSize),legend.text.align = 0)
p1 = p1 + scale_colour_manual(values = c("dodger blue","deep sky blue"),labels = c(expression(g == 6.3),expression(g == 13.4)))

p1 = p1 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps3f = p1
# ggsave(filename = file.path(figpath,sprintf("SuppFig3F.pdf",fstem_oof)),plot=p1)
p1

# Put plots together
p_final = (ps3a | ps3b | ps3c) / (ps3d | ps3e | ps3f)
ggsave(filename = file.path(figpath,sprintf("SuppFig3.pdf",fstem_oof)), plot=p_final)
p_final

Supp Figure 4

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

# Load data
# H values
BOLD_H_g_rate_1_Magri_kernel_No_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_1.5_Magri_kernel_No_HPF.csv"))
BOLD_H_g_rate_2_Magri_kernel_No_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_2.0_Magri_kernel_No_HPF.csv"))
BOLD_H_g_rate_1_Can_kernel_No_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_1.5_Can_kernel_No_HPF.csv"))
BOLD_H_g_rate_2_Can_kernel_No_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_2.0_Can_kernel_No_HPF.csv"))
BOLD_H_g_rate_1_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_1.5_Can_kernel_Yes_HPF.csv"))
BOLD_H_g_rate_2_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"BOLD_H_g_rate_2.0_Can_kernel_Yes_HPF.csv"))

# Correlation of LFP power and BOLD
corr_matrix_Magri_kernel_No_HPF = read.csv(file.path(datapath,"corr_matrix_Magri_kernel_No_HPF.csv"))
alpha_band_Magri_kernel_No_HPF = read.csv(file.path(datapath,"alpha_band_Magri_kernel_No_HPF.csv"))
beta_band_Magri_kernel_No_HPF = read.csv(file.path(datapath,"beta_band_Magri_kernel_No_HPF.csv"))
gamma_band_Magri_kernel_No_HPF = read.csv(file.path(datapath,"gamma_band_Magri_kernel_No_HPF.csv"))
LFP_band_Magri_kernel_No_HPF = read.csv(file.path(datapath,"LFP_band_Magri_kernel_No_HPF.csv"))

corr_matrix_Can_kernel_No_HPF = read.csv(file.path(datapath,"corr_matrix_Can_kernel_No_HPF.csv"))
alpha_band_Can_kernel_No_HPF = read.csv(file.path(datapath,"alpha_band_Can_kernel_No_HPF.csv"))
beta_band_Can_kernel_No_HPF = read.csv(file.path(datapath,"beta_band_Can_kernel_No_HPF.csv"))
gamma_band_Can_kernel_No_HPF = read.csv(file.path(datapath,"gamma_band_Can_kernel_No_HPF.csv"))
LFP_band_Can_kernel_No_HPF = read.csv(file.path(datapath,"LFP_band_Can_kernel_No_HPF.csv"))

corr_matrix_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"corr_matrix_Can_kernel_Yes_HPF.csv"))
alpha_band_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"alpha_band_Can_kernel_Yes_HPF.csv"))
beta_band_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"beta_band_Can_kernel_Yes_HPF.csv"))
gamma_band_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"gamma_band_Can_kernel_Yes_HPF.csv"))
LFP_band_Can_kernel_Yes_HPF = read.csv(file.path(datapath,"LFP_band_Can_kernel_Yes_HPF.csv"))

cLims = c(-0.06,0.28)
cMid = 0.125
# cLims = round(range(c(corr_matrix_Magri_kernel_No_HPF$r,
#                 corr_matrix_Can_kernel_No_HPF$r,
#                 corr_matrix_Can_kernel_Yes_HPF$r)),2)
# Supp Fig 4A -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

colnames(corr_matrix_Magri_kernel_No_HPF) <- c('X.1','X','Y','r')

# jet colormap
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

p = ggplot(corr_matrix_Magri_kernel_No_HPF, aes(X, Y, fill= r)) + geom_tile()
# p = p + scale_fill_gradientn(colors = jet.colors(7),breaks = seq(0,0.25,0.05))
# p = p + scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), 
                             # breaks = seq(0,0.25,0.05))
p = p + scale_fill_gradient2(low="blue", mid="white", high="red",
                             midpoint=cMid, limits=cLims)
p = p + xlab("Delay (s)") + ylab("Frequency (Hz)")+ ggtitle("Correlation between LFP power and BOLD")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + scale_y_reverse(breaks = round(seq(100, 0, by = -20),1))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.title = element_text(size=fontSize,hjust=0.2),
        legend.key.size = unit(3, "mm"))

ps4a = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4A.pdf",fstem_oof)))
p

# Supp Fig 4B -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
alpha_band_Magri_kernel_No_HPF["variable"] <- 'alpha'
beta_band_Magri_kernel_No_HPF["variable"] <- 'beta'
gamma_band_Magri_kernel_No_HPF["variable"] <- 'gamma'
LFP_band_Magri_kernel_No_HPF["variable"] <- 'LFP'

merged_data <- rbind(alpha_band_Magri_kernel_No_HPF,
                     beta_band_Magri_kernel_No_HPF,
                     gamma_band_Magri_kernel_No_HPF,
                     LFP_band_Magri_kernel_No_HPF)

p = ggplot()
p = p + geom_line(data = merged_data, aes(x = time, y = band, colour = variable),size = 0.5)
p = p + xlab("Delay (s)") + ylab("Correlation") + ggtitle("Correlation in selected frequency bands")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.9, 0.85),legend.text = element_text(size=fontSize),legend.text.align = 0)
p = p + scale_colour_manual(values = c("#006400","red","blue","black"))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps4b = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4B.pdf",fstem_oof)))
p

# Supp Fig 4C -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

# Box plot
box1 <- data.frame(seq(1,60,1),seq(1,60,1))
box2 <- data.frame(seq(1,60,1),seq(1,60,1))

# Clustering of H in 3 groups
box1[1:20,2] = BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(0,60,3),3]
box1[21:40,2] = BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(1,60,3),3]
box1[41:60,2] = BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(2,60,3),3]
box1[1:20,1] = format(round(mean(BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box1[21:40,1] = format(round(mean(BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box1[41:60,1] = format(round(mean(BOLD_H_g_rate_1_Magri_kernel_No_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box1) <- c("var","metric")

box2[1:20,2] = BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(0,60,3),3]
box2[21:40,2] = BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(1,60,3),3]
box2[41:60,2] = BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(2,60,3),3]
box2[1:20,1] = format(round(mean(BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box2[21:40,1] = format(round(mean(BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box2[41:60,1] = format(round(mean(BOLD_H_g_rate_2_Magri_kernel_No_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box2) <- c("var","metric")

# boxplot
p2 <-ggplot(box1, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p2 = p2 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p2 = p2 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p2 = p2 + annotate(geom="text", x=3, y=1.8, label="1.5 sp/s", color="red", size = 6)

# Other props.
p2 = p2 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.85) + ggtitle("1.5 sp/s")
p2 = p2 + theme(plot.title = element_text(hjust = 0.5))
#p2 = p2 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p2 = p2 + scale_color_manual(values=c("#808000", "#808000", "#808000"))
p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

# boxplot
p3 <-ggplot(box2, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p3 = p3 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p3 = p3 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p3 = p3 + annotate(geom="text", x=3, y=1.8, label="2 sp/s", color="red", size = 6)

# Other props.
p3 = p3 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.85)+ ggtitle("2 sp/s")
p3 = p3 + theme(plot.title = element_text(hjust = 0.5))
#p3 = p3 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p3 = p3 + scale_color_manual(values=c("#FA8072", "#FA8072", "#FA8072"))
p3 = p3 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps4c = (p2 | p3)
# ggsave(filename = file.path(figpath,sprintf("SuppFig4C.pdf",fstem_oof)),plot=ps4c)
ps4c

# Supp Fig 4D -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

colnames(corr_matrix_Can_kernel_No_HPF) <- c('X.1','X','Y','r')

# jet colormap
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

p = ggplot(corr_matrix_Can_kernel_No_HPF, aes(X, Y, fill= r)) + geom_tile()
# p = p + scale_fill_gradientn(colors = jet.colors(7),breaks = seq(0,0.25,0.05))
# p = p + scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100),
#                              breaks = seq(0,0.25,0.05))
p = p + scale_fill_gradient2(low="blue", mid="white", high="red", 
                       midpoint=cMid, limits=cLims)
p = p + xlab("Delay (s)") + ylab("Frequency (Hz)")+ ggtitle("Correlation between LFP power and BOLD")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + scale_y_reverse(breaks = round(seq(100, 0, by = -20),1))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.title = element_text(size=fontSize,hjust=0.2),
        legend.key.size = unit(3, "mm"))

ps4d = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4D.pdf",fstem_oof)))
p

#Supp Fig 4E ------------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
alpha_band_Can_kernel_No_HPF["variable"] <- 'alpha'
beta_band_Can_kernel_No_HPF["variable"] <- 'beta'
gamma_band_Can_kernel_No_HPF["variable"] <- 'gamma'
LFP_band_Can_kernel_No_HPF["variable"] <- 'LFP'

merged_data <- rbind(alpha_band_Can_kernel_No_HPF,
                     beta_band_Can_kernel_No_HPF,
                     gamma_band_Can_kernel_No_HPF,
                     LFP_band_Can_kernel_No_HPF)

p = ggplot()
p = p + geom_line(data = merged_data, aes(x = time, y = band, colour = variable),size = 0.5)
p = p + xlab("Delay (s)") + ylab("Correlation") + ggtitle("Correlation in selected frequency bands")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
p = p + theme(legend.title = element_blank(),legend.position = c(0.9, 0.85),legend.text = element_text(size=fontSize),legend.text.align = 0)
p = p + scale_colour_manual(values = c("#006400","red","blue","black"))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps4e = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4E.pdf",fstem_oof)))
p

# Supp Fig 4F -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

# Box plot
box1 <- data.frame(seq(1,60,1),seq(1,60,1))
box2 <- data.frame(seq(1,60,1),seq(1,60,1))

# Clustering of H in 3 groups
box1[1:20,2] = BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(0,60,3),3]
box1[21:40,2] = BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(1,60,3),3]
box1[41:60,2] = BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(2,60,3),3]
box1[1:20,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box1[21:40,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box1[41:60,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_No_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box1) <- c("var","metric")

box2[1:20,2] = BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(0,60,3),3]
box2[21:40,2] = BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(1,60,3),3]
box2[41:60,2] = BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(2,60,3),3]
box2[1:20,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box2[21:40,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box2[41:60,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_No_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box2) <- c("var","metric")

# boxplot
p2 <-ggplot(box1, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p2 = p2 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p2 = p2 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p2 = p2 + annotate(geom="text", x=3, y=1.9, label="1.5 sp/s", color="red", size = 6)

# Other props.
p2 = p2 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.95) + ggtitle("1.5 sp/s")
p2 = p2 + theme(plot.title = element_text(hjust = 0.5))
#p2 = p2 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p2 = p2 + scale_color_manual(values=c("#808000", "#808000", "#808000"))
p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

# boxplot
p3 <-ggplot(box2, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p3 = p3 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p3 = p3 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p3 = p3 + annotate(geom="text", x=3, y=1.9, label="2 sp/s", color="red", size = 6)

# Other props.
p3 = p3 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.95) + ggtitle("2 sp/s")
p3 = p3 + theme(plot.title = element_text(hjust = 0.5))
#p3 = p3 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p3 = p3 + scale_color_manual(values=c("#FA8072", "#FA8072", "#FA8072"))
p3 = p3 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps4f = (p2 | p3)
# ggsave(filename = file.path(figpath,sprintf("SuppFig4F.pdf",fstem_oof)),plot=ps4f)
ps4f

# Supp Fig 4G -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

colnames(corr_matrix_Can_kernel_Yes_HPF) <- c('X.1','X','Y','r')

# jet colormap
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

p = ggplot(corr_matrix_Can_kernel_Yes_HPF, aes(X, Y, fill= r)) + geom_tile()
# p = p + scale_fill_gradientn(colors = jet.colors(7), breaks = seq(0,0.25,0.05))
# p = p + scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100),
#                              breaks = seq(0,0.25,0.05))
p = p + scale_fill_gradient2(low="blue", mid="white", high="red", 
                       midpoint=cMid, limits=cLims)
p = p + xlab("Delay (s)") + ylab("Frequency (Hz)")+ ggtitle("Correlation between LFP power and BOLD")
p = p + theme(plot.title = element_text(hjust = 0.5))
p = p + scale_y_reverse(breaks = round(seq(100, 0, by = -20),1))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.title = element_text(size=fontSize,hjust=0.2),
        legend.key.size = unit(3, "mm"))

ps4g = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4G.pdf",fstem_oof)))
p

# Supp Fig 4H -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
alpha_band_Can_kernel_Yes_HPF["variable"] <- 'alpha'
beta_band_Can_kernel_Yes_HPF["variable"] <- 'beta'
gamma_band_Can_kernel_Yes_HPF["variable"] <- 'gamma'
LFP_band_Can_kernel_Yes_HPF["variable"] <- 'LFP'

merged_data <- rbind(alpha_band_Can_kernel_Yes_HPF,
                     beta_band_Can_kernel_Yes_HPF,
                     gamma_band_Can_kernel_Yes_HPF,
                     LFP_band_Can_kernel_Yes_HPF)

p = ggplot()
p = p + geom_line(data = merged_data, aes(x = time, y = band, colour = variable),size = 0.5)
p = p + xlab("Delay (s)") + ylab("Correlation") + ggtitle("Correlation in selected frequency bands")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.9, 0.85),legend.text = element_text(size=fontSize),legend.text.align = 0)
p = p + scale_colour_manual(values = c("#006400","red","blue","black"))

p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps4h = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig4H.pdf",fstem_oof)))
p

# Supp Fig 4I -----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)

# Box plot
box1 <- data.frame(seq(1,60,1),seq(1,60,1))
box2 <- data.frame(seq(1,60,1),seq(1,60,1))

# Clustering of H in 3 groups
box1[1:20,2] = BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(0,60,3),3]
box1[21:40,2] = BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(1,60,3),3]
box1[41:60,2] = BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(2,60,3),3]
box1[1:20,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box1[21:40,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box1[41:60,1] = format(round(mean(BOLD_H_g_rate_1_Can_kernel_Yes_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box1) <- c("var","metric")

box2[1:20,2] = BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(0,60,3),3]
box2[21:40,2] = BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(1,60,3),3]
box2[41:60,2] = BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(2,60,3),3]
box2[1:20,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(0,60,3),2]), 2), nsmall = 2)
box2[21:40,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(1,60,3),2]), 2), nsmall = 2)
box2[41:60,1] = format(round(mean(BOLD_H_g_rate_2_Can_kernel_Yes_HPF[seq(2,60,3),2]), 2), nsmall = 2)
colnames(box2) <- c("var","metric")

# boxplot
p2 <-ggplot(box1, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p2 = p2 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p2 = p2 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p2 = p2 + annotate(geom="text", x=3, y=1.9, label="1.5 sp/s", color="red", size = 6)

# Other props.
p2 = p2 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.85) + ggtitle("1.5 sp/s")
p2 = p2 + theme(plot.title = element_text(hjust = 0.5))
#p2 = p2 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p2 = p2 + scale_color_manual(values=c("#808000", "#808000", "#808000"))
p2 = p2 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

# boxplot
p3 <-ggplot(box2, aes(x = fct_reorder(var,metric,.desc=TRUE), y = metric, colour = var)) + 
  geom_boxplot(fill = NA, outlier.shape = NA)
# Box plot with dot plot
p3 = p3 + geom_jitter(shape=16, position=position_jitter(0.2))
# #  Add p-value
# p3 = p3 + stat_compare_means(method = "anova",size = 5)
# # Annotate external rate
# p3 = p3 + annotate(geom="text", x=3, y=1.9, label="2 sp/s", color="red", size = 6)

# Other props.
p3 = p3 + xlab( expression(paste( g == g[I]/g[E])) ) + ylab("H") + ylim(1.5,1.85) + ggtitle("2 sp/s")
p3 = p3 + theme(plot.title = element_text(hjust = 0.5))
#p3 = p3 + scale_y_discrete(breaks = round(seq(1.4, 1.7, by = 0.1),1))
p3 = p3 + scale_color_manual(values=c("#FA8072", "#FA8072", "#FA8072"))
p3 = p3 + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps4i = (p2 | p3)
# ggsave(filename = file.path(figpath,sprintf("SuppFig4I.pdf",fstem_oof)),plot=ps4f)
ps4i

# Put plots together
p_final = ((ps4a | ps4b | ps4c) + plot_layout(widths=c(1,1,2.25))) / 
  ((ps4d | ps4e | ps4f) + plot_layout(widths=c(1,1,2.25))) / 
  ((ps4g | ps4h | ps4i) + plot_layout(widths=c(1,1,2.25))) 
ggsave(filename = file.path(figpath,sprintf("SuppFig4.pdf",fstem_oof)), plot=p_final)
p_final

Supp Figure 5

# Global properties of plots
fontSize = 8
fstem_oof = "oof"

legend_labels_slope = c(expression(paste(nu[0] == 1.5," sp/s")),
                        expression(paste(nu[0] == 2," sp/s")),
                        'Baseline')

# Load data
# Vth
# LFP slopes 
lf_slopes_vth_rate_1 = read.csv(file.path(datapath,"lf_slopes_vth_rate_1.5.csv"))
lf_slopes_vth_rate_2 = read.csv(file.path(datapath,"lf_slopes_vth_rate_2.0.csv"))
hf_slopes_vth_rate_1 = read.csv(file.path(datapath,"hf_slopes_vth_rate_1.5.csv"))
hf_slopes_vth_rate_2 = read.csv(file.path(datapath,"hf_slopes_vth_rate_2.0.csv"))

# EL
# LFP slopes 
lf_slopes_EL_rate_1 = read.csv(file.path(datapath,"lf_slopes_EL_rate_1.5.csv"))
lf_slopes_EL_rate_2 = read.csv(file.path(datapath,"lf_slopes_EL_rate_2.0.csv"))
hf_slopes_EL_rate_1 = read.csv(file.path(datapath,"hf_slopes_EL_rate_1.5.csv"))
hf_slopes_EL_rate_2 = read.csv(file.path(datapath,"hf_slopes_EL_rate_2.0.csv"))

# Vth
fooof_ap_vth_rate_1 = read.csv(file.path(datapath,"fooof_ap_vth_rate_1.5.csv"))
fooof_ap_vth_rate_2 = read.csv(file.path(datapath,"fooof_ap_vth_rate_2.0.csv"))

# EL
fooof_ap_EL_rate_1 = read.csv(file.path(datapath,"fooof_ap_EL_rate_1.5.csv"))
fooof_ap_EL_rate_2 = read.csv(file.path(datapath,"fooof_ap_EL_rate_2.0.csv"))

# Baseline lines
baseline_vth_slopes_1 <- data.frame("var" = -52,"metric"=seq(-1, 1, length.out=100),"variable"='Baseline')
baseline_vth_slopes_2 <- data.frame("var" = -52,"metric"=seq(-2.75, -1, length.out=100),"variable"='Baseline')

baseline_EL_slopes_1 <- data.frame("var" = -70,"metric"=seq(-1, 1, length.out=100),"variable"='Baseline')
baseline_EL_slopes_2 <- data.frame("var" = -70,"metric"=seq(-2.75, -1, length.out=100),"variable"='Baseline')

baseline_vth_slopes <- data.frame("var" = -52,"metric"=seq(-1.25, 0.5, length.out=100),"variable"='Baseline')
baseline_EL_slopes <- data.frame("var" = -70,"metric"=seq(-1.25, 0.5, length.out=100),"variable"='Baseline')

# Supp Fig 5A1 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
lf_slopes_vth_rate_1["variable"] <- '1.5 spikes/s'
lf_slopes_vth_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(lf_slopes_vth_rate_1,lf_slopes_vth_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_vth_slopes_1, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_vth_slopes_1, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( V[th]," (mV)")) ) + ylab("1/f slope") + ggtitle("Low-frequencies <30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072","black"),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps5a1 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5A1.pdf",fstem_oof)))
p

# Supp Fig 5A2 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
hf_slopes_vth_rate_1["variable"] <- '1.5 spikes/s'
hf_slopes_vth_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(hf_slopes_vth_rate_1,hf_slopes_vth_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_vth_slopes_2, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_vth_slopes_2, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste(  V[th]," (mV)")) ) + ylab("1/f slope") + ggtitle("High-frequencies >30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072","black"),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps5a2 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5A2.pdf",fstem_oof)))
p

# Supp Fig 5A3 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
fooof_ap_vth_rate_1["variable"] <- '1.5 spikes/s'
fooof_ap_vth_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(fooof_ap_vth_rate_1,fooof_ap_vth_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_vth_slopes, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_vth_slopes, aes(x = var, y = metric, colour = FALSE))
p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( V[th]," (mV)")) ) + ylab("1/f slope") + ggtitle("Aperiodic slope")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize))

ps5a3 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5A3.pdf",fstem_oof)))
p

# Supp Fig 5B1 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
lf_slopes_EL_rate_1["variable"] <- '1.5 spikes/s'
lf_slopes_EL_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(lf_slopes_EL_rate_1,lf_slopes_EL_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_EL_slopes_1, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_EL_slopes_1, aes(x = var, y = metric, colour = FALSE))
# p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( E[L]," (mV)")) ) + ylab("1/f slope") + ggtitle("Low-frequencies <30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072","black"),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps5b1 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5B1.pdf",fstem_oof)))
p

# Supp Fig 5B2 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
hf_slopes_EL_rate_1["variable"] <- '1.5 spikes/s'
hf_slopes_EL_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(hf_slopes_EL_rate_1,hf_slopes_EL_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_EL_slopes_2, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_EL_slopes_2, aes(x = var, y = metric, colour = FALSE))
# p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste(  E[L]," (mV)")) ) + ylab("1/f slope") + ggtitle("High-frequencies >30 Hz")
p = p + theme(plot.title = element_text(hjust = 0.5))
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072",'black'),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = "none")

ps5b2 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5B2.pdf",fstem_oof)))
p

# Supp Fig 5B3 ----------------------------------------------------------------
labelidx2use = c(1,6,11,16,21)
fooof_ap_EL_rate_1["variable"] <- '1.5 spikes/s'
fooof_ap_EL_rate_2["variable"] <- '2 spikes/s'
merged_data <- rbind(fooof_ap_EL_rate_1,fooof_ap_EL_rate_2)

p = ggplot()
p = p + geom_point(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        geom_smooth(data = merged_data, aes(x = var, y = metric, colour = variable)) +
        guides(data = merged_data, aes(x = var, y = metric, colour = FALSE))
p = p + geom_line(data = baseline_EL_slopes, aes(x = var, y = metric, colour = variable),linetype = "dotdash") +
        guides(data = baseline_EL_slopes, aes(x = var, y = metric, colour = FALSE))
#p = p + scale_x_reverse(breaks = round(seq(-51.5, -53, by = -0.5),1))
p = p + xlab( expression(paste( E[L]," (mV)")) ) + ylab("1/f slope") + ggtitle("Aperiodic slope")
p = p + theme(plot.title = element_text(hjust = 0.5)) + guides(colour=FALSE)
# p = p + theme(legend.title = element_blank(),legend.position = c(0.8, 0.2),legend.text = element_text(size=fontSize),legend.text.align = 0)
# p = p + guides(color = guide_legend(override.aes = list(shape= NA,linetype = c("solid","solid","dotdash"),size=c(1,1,1))))
p = p + scale_colour_manual(values = c("#808000","#FA8072","black"),labels = legend_labels_slope)
p = p + theme(axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text.x = element_text(size=fontSize),
        plot.title = element_text(size=fontSize),
        legend.position = 'none')

ps5b3 = p
# ggsave(filename = file.path(figpath,sprintf("SuppFig5B3.pdf",fstem_oof)))
p

# Put plots together
p_final = (ps5a1 | ps5a2 | ps5a3) / (ps5b1 | ps5b2 | ps5b3)
ggsave(filename = file.path(figpath,sprintf("SuppFig5.pdf",fstem_oof)), plot=p_final)
p_final